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ABSTRACT 

We investigated the effect of inverse Compton scattering in mildly relativistic 
static and moving plasmas with low optical depth using Monte Carlo simulations, and 
calculated the Sunyaev-ZePdovich effect in the cosmic background radiation. Our 
semi- analytic method is based on a separation of photon diffusion in frequency and 
real space. We use Monte Carlo simulation to derive the intensity and frequency 
of the scattered photons for a monochromatic incoming radiation. The outgoing 
spectrum is determined by integrating over the spectrum of the incoming radiation 
using the intensity to determine the correct weight. This method makes it possible 
to study the emerging radiation as a function of frequency and direction. As a first 
application we have studied the effects of finite optical depth and gas infall on the 
Sunyaev-ZePdovich effect (not possible with the extended Kompaneets equation) and 
discuss the parameter range in which the Boltzmann equation and its expansions 
can be used. For high temperature clusters (ksT e ^ 15 keV) relativistic corrections 
based on a fifth order expansion of the extended Kompaneets equation seriously 
underestimate the Sunyaev-ZePdovich effect at high frequencies. The contribution 
from plasma infall is less important for reasonable velocities. We give a convenient 
analytical expression for the dependence of the cross-over frequency on temperature, 
optical depth, and gas infall speed. Optical depth effects are often more important 
than relativistic corrections, and should be taken into account for high-precision work, 
but are smaller than the typical kinematic effect from cluster radial velocities. 
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1. Introduction 

Inverse Compton scattering of the cosmic microwave background radiation (CMBR) by hot 
electrons in the atmospheres of clusters of galaxies, the Sunyaev-Zel’dovich (SZ) effect (Sunyaev 
and Zel’dovich 1980), has become a powerful tool in astrophysics. It is one of the most important 
secondary effects which cause fluctuations in the CMBR. We will refer to the effect arising 
from static gas as the static SZ (SSZ) effect, and that arising from gas with bulk motion as the 
kinematic SZ (KSZ) effect. Fluctuations in the CMBR caused by the SZ effects in an ensemble of 
clusters of galaxies should dominate on angular scales less than few arcminutes. The nature of 
these fluctuations depends on the evolution of clusters, and so it is a test of structure formation 
theories (eg. Aghanim et al. 1998; Molnar and Birkinshaw 1998). 

Observations of the SSZ effect, begun in the 1970s, have now become routine in the 90s 
with dedicated instruments using the latest receiver technology (for reviews see Rephaeli 1995b; 
Birkinshaw 1998). SZ effect and x-ray measurements probe the physical conditions in the 
intracluster gas in clusters of galaxies, and allow us to deduce the distance to the cluster without 
additional assumptions. This provides a useful independent method for determining the Hubble 
constant. 

Observations of the KSZ effect are much more difficult since it is typically an order of 
magnitude smaller than the SSZ effect and has the same spectrum as the primordial fluctuations 
in the CMBR. The KSZ effect provides a method of measuring radial peculiar velocities of 
clusters, and even without the tangential velocity component, which might be determined using 
the Rees-Sciama (RS) effect (Rees and Sciama 1968; Birkinshaw and Gull 1983; Gurvits and 
Mitrofanov 1986; Aghanim et al. 1998; Molnar and Birkinshaw 1998) it should provide important 
information on large scale velocity fields, which are closely related to the large scale density 
distributions and thus to the average total mass density in the Universe. Useful limits on the size 
of the KSZ effect for two clusters have recently been reported by Holzapfel et al. (1997a). 

Most discussions of the SSZ effect have been based on non-relativistic calculations of its 
amplitude, made via a Fokker-Planck type expansion of the Boltzmann equation (Kompaneets 
1957). The advantage of this approach is that in interesting cases (electron temperature, T e , 
much greater than the temperature of the incoming radiation) it provides a convenient analytical 
solution for the spectrum of the emerging radiation. However, it has been recognized recently 
that relativistic effects become important for clusters with kT e ^ 10 keV. Rephaeli (1995a) 
provided a relativistic solution for the SSZ effect as a series expansion in the optical depth (*C 1 
for clusters). There is no exact analytical solution. The numerical integrals involved axe tractable 
in the single-scattering approximation, which is usually adequate in clusters. Faxgion, Konoplich 
and Salis (1996) developed exact expressions for relativistic inverse Compton scattering of a laser 
beam with monochromatic isotropic radiation, finding good agreement with the approximations 
of Jones (1968). Their expression for the frequency redistr bution function (FRDF) for scattering 
of monoenergetic electrons with monochromatic photons agrees with Rephaeli’s result. Rephaeli’s 
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method has been used in a series of papers to evaluate the SSZ effect for hot clusters. Relativistic 
corrections to the SSZ effect, and to the standard thermal bremsstrahlung formulae, were applied 
to determine the Hubble constant in hot clusters by Rephaeli and Yankovich (1997), however, 
their corrections of the thermal bremsstrahlung equation were further corrected by Hughes and 
Birkinshaw (1998). Holzapfel et al. (1997b) used the relativistic results in their determination of 
the Hubble constant from observations of cluster Abell 2163. 

The most general treatment of Compton scattering in static and moving media has been 
derived by Psaltis and Lamb (1997) as a series expansion. As Challinor and Lasenby (1998b) 
noted, however, more terms in the expansion should be taken into account for accurate treatment 
of clusters of galaxies. Recently the Kompaneets equation has been extended to contain relativistic 
corrections to the SSZ and KSZ effects (Stebbins 1997; Challinor and Lasenby 1998a, b; Itoh, 
Kohyama and Nozawa 1998; Nozawa, Itoh and Kohyama 1998; Sazonov and Sunyaev 1998b). 
Starting from the Boltzmann equation, an expansion in the small parameters of the dimensionless 
temperature, 0 € = fcT e /(m e c 2 ), fractional energy change in a scattering, {hv* — hv)jksT e ^ and 
dimensionless radial velocity for the KSZ effect, u ra( */c, leads to a Fokker-Planck type equation 
(the extended Kompaneets equation). Corrections up to the fifth order in 0 e have been derived 
(Itoh et al. 1998). These calculations demonstrate the importance of the relativistic effects (in 
accordance with the results of Rephaeli 1995a). Note however, that, as Challinor and Lasenby 
(1998a) emphasized, the extended Kompaneets equation is a result of an asymptotic series 
expansion, therefore it is important to estimate the validity of the expansion using other methods. 
Nozawa et al. compared the convergence of their expansion to a direct numerical evaluation of 
the Boltzmann collision integral, and concluded that in the Rayleigh- Jeans region the relativistic 
corrections give accurate results in the entire range of cluster temperatures. Significant deviations 
axe found at higher frequencies for high temperature clusters. 

The SSZ and KSZ effects must be separated in order to extract information on peculiar 
velocities. Fortunately the two effects have different frequency dependence. The maximum of 
the KSZ effect (in thermodynamic temperature units) occurs at about the “cross-over” frequency 
where the SSZ effect changes sign from being a decrement to an increment. In a non-relativistic 
treatment the cross-over frequency is a constant, 218 GHz, independent of electron temperature, 
optical depth, and all other parameters. Rephaeli (1995a) showed that in the relativistic case 
the cross-over frequency depends on the temperature, and his results were used by Holzapfel 
et al. (1997a) in determining peculiar velocities of two clusters. Sazonov and Sunyaev (1998b) 
and Nozawa et al. (1998) give approximations for the cross-over frequency as a function of 
dimensionless temperature and radial peculiar velocity. They also conclude that relativistic 
corrections to the cross-over frequency are important, and should be taken into account in future 
experiments. 

Other methods have been used to investigate inverse Compton scattering, such as numerical 
integration of the collision integral (Corman 1970), multiple scattering methods (Wright 1979), 
and Monte Carlo simulations. Simulations of inverse Compton scattering in relativistic and 
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non-relativistic plasma have been carried out for embedded sources (Pozdnyakov, Sobol, and 
Sunyaev 1983; Haardt and Maraschi 1993; Hua and Titarchuk 1995). Gull and Garret (1998) 
used Monte Carlo methods to evaluate the Boltzmann colllsional integral. Sazonov and Sunyaev 
(1998a) used Monte Carlo simulations to derive the SZ thermal and kinematic effects. 

In this paper we study the effect of optical depth and non-uniform bulk motion on the SZ 
effect using a Monte Carlo method to calculate the frequency redistribution function. The inverse 
Compton scattering of CMBR photons is treated in the Thomson limit for static and infalling 
plasmas (SSZ and KSZ effects) with spherical symmetry, uniform density distribution, and low 
optical depth over a wide range of gas temperatures and observed frequency. We apply our results 
to clusters of galaxies assuming a static and radially infalling (or collapsing) gas component. 


2. The Method 
2.1. Formsdism 

The emerging intensity of a beam of radiation in the line of sight after passage through a 
scattering atmosphere can be expressed as 


I{x) = e T B(x) + I ac {x), (1) 

where the incoming intensity (hereafter assumed to be Planckian), B(x), is attenuated by a 
factor depending on the optical depth r and we add the scattered intensity into the beam, I sc , at 
dimensionless frequency x = hv/ksTcB, where h, u, ks and Tcb are the Planck constant, the 
frequency, the Boltzmann constant, and the temperature of the CMBR, T C b — 2.728 ± 0.002 K 
(Fixsen et al. 1996). The scattered intensity may be determined by using 

Isc{v) = J P{a) B ac {i/ 0 )ds, ( 2 ) 

where the intensity of the scattered incoming radiation B ac ( v) = w ac B(i /), and w ac is the fraction 
of photons scattered into the beam, the weight for this beam. P(s) is the frequency redistribution 
function (FRDF), which specifies the probability of scattering from uq to v as a function of 
the logarithm of the dimensionless frequency s = \d.(v/vq). This function gives the frequency 
distribution (normalized to unity) of the scattered radiation (we give P(s) a slightly different 
definition from the one given by Birkinshaw 1998). This decomposition is possible because in our 
approximation the fractional frequency change is independent of the frequency (see equation 16 
later). The change of the intensity in the line of sight may be expressed as 


AI(x) = (e~ T -l)B(x) + I ac (x). 


( 3 ) 
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For isotropic and homogeneous scattering conditions, so that the scattering parameters do not 
depend on where the scattering happens, the weight 


w sc = (1 - e T ), (4) 

which means that the out-scattered radiation is balanced by the same amount of in-scattered 
radiation, and so there would be no net intensity change (A I(x) = 0) if there were no frequency 
change ( P(s ) is the Dirac delta function). This assumption breaks down, for example, when 
there is relativistic bulk motion, which introduces anisotropy in the scattering via the relativistic 
beaming effect. Our task is to determine P(s) and w sc . We used equation (4) as a check for our 
Monte Carlo simulations for static models (where it must be correct), and verified that it is a good 
approximation for our models with infall by counting outgoing photons. We derive P(s) using 
a Monte Carlo method. At low optical depth (r ^ 1), the problem is suitable for Monte Carlo 
simulation because we do not have to follow photons through many scatterings in the medium (the 
average number of scatterings being r). Although in this limit most of the photons do not scatter, 
and hence provide no information on P(s), this is not a problem since we can use the method of 
forced first scattering (see below). 


2.2. Monte Carlo Method 

We give a short description of the method here, for a more detailed description see Molnar 
(1998). 

We assume an isotropic incoming low temperature radiation field (the CMBR). We use forced 
first scatterings to study the inverse Compton process. The photons Compton scatter from an 
electron population with a relativistic Maxwellian distribution of momenta in the rest frame of 
bulk motion. We compute scattering probabilities in the rest frame of the electron in the Thomson 
limit (Chandrasekhar 1950). This involves coordinate transformations from the observer’s frame 
to the rest frame of the bulk motion and to the rest frame of the electron. We assume time 
translation invariance and spherical symmetry. Time translational invariance is not exact for 
our model with infall, so that we make a snap-shot approximation. The error arising from this 
approximation is less than the light crossing time over the infall time (« v r /c = /3 r ), which is only 
a few per cent of the infall term for our models. 

In most cases we use the inverse method to generate the desired probability distribution. We 
use a rejection method when the inverse method leads to non-invertible functions or is too slow: 
for a general description of generating probability distributions cf. Pozdnyakov et al. (1983); Press 
et al. (1992). The former reference describes an alternative Monte Carlo method to treat inverse 
Compton scattering. 

In the description that follows we use the word “photon” in the singular to refer to one Monte 
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Carlo “photon”, one experiment in our simulation. We use a weight, ttfi n , to express the number 
of photons this one experiment represents (the weight does not have to be an integer). We carried 
out the simulation in five steps. 

Step 1. The position and direction of incoming photons: 

We assume that the photons arrive uniformly on a unit sphere (the radius of the gas is scaled to 
unity). In a coordinate system which is placed at the point of impact, the direction cosine of the 
incoming photons from the normal, can be sampled as 

Min = (5) 

where we use [RN] to indicate a uniformly distributed random number drawn each time when it 
occurs. The azimuthal angle is assumed to be uniformly distributed between 0 and 2i r. 

Step 2 . distance to the forced first scattering : 

We use forced first scattering, which means that we take the probability of scattering equal to one 
on the photons’ original line of flight through the gas. Using the inversion method, the optical 
depth at which the incoming photon scatters becomes 


ri = - ln(l - [RN] • (1 - e~ T ™ ax )), (6) 

where r^ ax is the maximum optical depth in the line of sight (the superscript refers to the 
number of times the photon has already scattered; zero in this case). Since we are using 
forced first collisions, we have to account for the fraction of photons which are unscattered 
on their path through the cloud. The scattered weight may be obtained from equation (4): 
w sc = (1 — e“ T ™*)u;i n , and stays the same during subsequent (unforced) scatterings. The weight 
of the photons passing through the cloud without scattering is w^t — Wi n e~ T ™ ax , 

Step 3. scattering: 

At the calculated position of the (forced first) scattering, we use the direction of the incoming 
photon and obtain the direction of propagation and frequency of the scattered photon. In the 
case that the gas is moving, we make a Lorentz transformation into the rest frame of the moving 
plasma. We sample the scattered electron’s dimensionless velocity, f3 e = v c /c, from a relativistic 
Maxwellian distribution 


P(f3 e ) d(3 e = N rel 7 5 e“ 7/e ' dp e , (7) 

where the normalization is 


N rel = (e e tf 2 (i/e e )) \ 


(8) 
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K 2 is the second order modified Bessel function of the second kind, and the dimensionless electron 
temperature is 


©e 


We used the rejection method to sample /3 e . 


fcj? Te 
m e c 2 ' 


( 9 ) 


The distribution of the direction of electron momenta is simplest in a frame in which the 
photon momentum unit vector points into one of the coordinate axes, z for example. In this 
coordinate system the probability distribution of /x e , the cosine of the angle between the unit 
vector of the direction of photon propagation (z axis) and electron velocity, is 


P(l*e,Pe) = ( l-0elle)/2. 


The inversion method leads to sampling fi e as 


( 10 ) 


' ,<= s:( 1 -l/ 1_2A ( 2|iW1 - 1 'f))’ (11) 

where the sign of the square root was determined so that in the limit of small electron velocities we 
recover the result for an isotropic distribution (/u e = 2 [RN] — 1). The angular distribution of the 
scattered electrons in the plane perpendicular to the momentum vector of the photon is isotropic 
(at azimuthal angle uniformly distributed between 0 and 2i r). In the rest frame of the electron, 
the cosine of the polar angle of the incoming photon is derived from a Lorentz transformation as 


_ — + fle 

fi ~ 1-/W 


( 12 ) 


where the negative sign in front of fi e is appropriate for an incoming photon. In the electron’s rest 
frame the scattering probability of a photon coming in with direction cosine and leaving with 
direction cosine y! is given by Chandrasekhar (1950) 


— g(l + M 2 ^ /2 + 2 ^ ” / i2 )( 1 ~ t 1 ' 2 ))' (13) 

/i' can be sampled using a uniform probability distribution by inversion of 


[JW]- = + (14) 

which leads to a cubic equation for 


(m 2 - ly 3 + (3 - M 2 y — | (2 [RN] + |) = 0. 


(15) 
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This cubic equation has a single real solution with absolute value of /i' less or equal to one. 

We now transfer the direction and frequency of the scattered radiation back to the observer’s 
frame. The dimensionless outgoing frequency of the photon normalized to the incoming frequency 
(expressed with the s parameter) becomes 






Pei*ej H" Pel* • 


(16) 


Step 4 - loop over scatterings: 

Having the point of scattering, the scattered frequency, and the direction of the scattered photon, 
we now sample the optical depth to the next scattering. We do not use forced scattering, so the 
optical depth follows from the usual (inverse) method as 


T n = “ ln[J?JV], (17) 

where the lower index on r refers to the nth scattering (n > 1), and we used [RN] = 1 — 
which is correct for uniform probability distributions. If r is less than the maximum optical depth 
in the direction of the photon momentum after the previous scattering, the photon is taken 

to have scattered within the cloud, and we calculate the new scattering direction and frequency 
of the scattered photon as in step 3. If Th ^ ? maxi the photon escaped (scattered n — 1 times), 
and we register the impact parameter (with the weight of the photon) and the dimensionless 
frequency, s, of the escaped photon. At the end of the simulation we sum the weights in every 
impact parameter bin to check our assumption of homogeneous scattering (equation 4). 

Step 5. The frequency redistribution function: 

When the simulation finishes, we derive the FRDF by binning the frequencies of the 
out-coming photons. The discrete probability distribution then becomes 


P(s k ) = N a Af k . (18) 

Mk is the number of photons in the k th bin, for which s* — As/2 < s < s* + As/2, where s* is the 
center of the k t/l bin, As is the width of the bin, and N s is the normalization, N s = 1/(AsIVmc)- 
Nmc is the number of Monte Carlo photons. This P(sfc) is our sampled approximation to P(s), 
and we then fit an exponential of polynomials to P(sjt) to get a convenient expression for the 
FRDF. Our fit gives an approximation accurate to better than half a percent, except in the (small) 
extended tails, where P(s) is under-represented. However, these regions lie 3 orders of magnitude 
below the peak, and the error arising from the fit is negligible for our applications. 
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2.3. Testing the code 

2.3.1. Single-scattering approximation 

We tested our code by comparing our single scattering Monte Carlo results for the FRDF 
(pMC) £ 0 thogg derived from Rephaeli (1995a). Rephaeli’s single-scattering approximation can be 
written as 


A(S) “ 32 N rel i 0 7 6 7/0 '^ 4(/l +h + h) d/?e ’ 


(19) 


where 


fl — ^ S 0e (^2 — M2 — Ml + Ml) 

h = fe - 3(4e* + 1) + 3e a (x 2 + X!) + ( g2 - Xl ) 

\ x 2 Xl / 

/ 3 = 2e s (e s + 1)(3 — Pi) In — , (20) 

X i 

where x\ = 1 — (3 e Mi, x? = 1 — (3 e j.i 2 - Mi and M 2 are defined by 



f-1 

s < 0 

Pl = 


s > 0 


1 

V 

1 

s < 0 

M2 = 

< Pe 



U 

s > 0 


( 21 ) 


and we used equation (16) to eliminate p/. This result (equation [19]) agrees with Fargion et al. 
(1996). These expressions can be integrated numerically, except when s = 0. In that case (3q = 0, 
and direct numerical integration is not possible because of the diverging /3~ 4 term. For small 
Pe {Pe < b — 0.1 for example) we can expand the logarithm, and use this expansion as a good 
approximation. For s = 0, equation (19) becomes 


Pi(>=0)= 8 < 22 > 

where the integrand is 

fo = 7e- 7/0 '/?- 4 (V Pi + 6 Pe + (3 - Pl) . (23) 

A Maclaurin expansion of the logarithm to order /3g gives an adequate approximation 
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/: * «■ - c i + s - § ■ + (i - ?k ■ + (hk) «■ < 24 > 

with no remaining divergent terms for the first integral in equation (22). 

We derive Pj MC (s) from our Monte Carlo simulation by using the results of only the first 
(forced) scatterings. Figures la and lb show our Monte Carlo results, P\ fC , superimposed on Pi 
from equations (19) and (20) for dimensionless temperatures 0 e = 0.03 and 0.3. The agreement is 
excellent, confirming that our Monte Carlo code is successfully reproducing Pi(s). 


2.3.2. Testing the Numerical Integral for the Intensity Change 

We derive the intensity change from the FRDF using a numerical integral (equation [3]). The 
Kompaneets approximation leads to the following FRDF: 


Pk 

where the Compton y parameter is 



(a ~ 3y) 2 \ 
4 y r 


(25) 


y = J n e GT ®edt, (26) 

where n e is the electron number density as a function of length in the line of sight measured by 
l and © e is the dimensionless temperature (for a discussion see for example Birkinshaw 1998; 
Molnar 1998). In order to check our numerical method, we used the Kompaneets FRDF (25) in 
the numerical integral in equation (3), and compared the resulting intensity change to that of 
obtained by the analytic solution for the Kompaneets approximation 


A r _ iox 4 e x ( e x + l A \ 
AIk V (e* - l) 2 V*e* - f / 


(27) 


where to = 2 (A:bTcb) 3 /(/ic) 2 ). We concluded that our numerical method is accurate better than 

0 . 1 %. 


2.3.3. Relativistic Corrections to Kompaneets equation 

We also compare our results to those from the extend' id Kompaneets equation up to the 5th 
order in 0 e (Itoh et al. 1998). Itoh et al. expressed the intensity change as 


A7 = 


An iox 3 
n e x — 1 ’ 


(28) 
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and provided expressions for ^ (note, that their y parameter is actually r, the optical depth). On 
Figure 2 we plot A I/t from the Kompaneets approximation, from Itoh et al.’s expansion, and for 
our single scattering Monte Carlo results. From the figure we conclude that our single scattering 
Monte Carlo result agrees with that of Itoh et al. at low temperatures, 6 e < 0.03 (T e < 15 keV). 
Deviations from the Itoh et al.’s result are already appearing at 0 e = 0.03, and become more 
pronounced at higher temperatures and high frequencies, as we would expect. 

We conclude that our simulation method passes these two tests, and can now be used to 
calculate the effects of multiple scattering and bulk velocity on the SZ effect. 


3. Results 

We performed a number of simulations for uniform density spherical models which are either 
static or have radial infalls with constant gas velocity at all radii. These models cover a range of 
to, the optical depth of zero impact parameter, electron temperature, T e , and gas infall speed, /3 r . 
Simulations with monochromatic incoming radiation were used to determine the FRDF. 

In Figure 3 we show the FRDF for the static model with r 0 = 0.05 at seven different 
temperatures. At higher temperatures more photons scatter into higher energies, thus the FRDFs 
are broader, and have lower peaks (since they are normalized to unity). We show the effect of 
finite optical depth in Figure 4: higher optical depth leads to more scatterings, and therefore more 
photons scattered to higher energies. Even for an optical depth as large as r 0 = 0.1 the change in 
the function is relatively small. In Figure 5 we show the effect on the FRDF of gas infall. Larger 
infall velocities cause more up-scattering of the photons, and hence more spreading of the FRDF, 
but the most obvious change is that the sharp peak at s = 0 is smoothed out by the motion of the 
plasma. At lower temperatures bulk motion causes larger departures from the static FRDF since 
the infall speed is larger relative to the electron thermal velocity. 

We used these results for the FRDF to calculate the spectra of scattered radiation by the 
convolution in equation (2). We derived the intensity change using equation (3). In the static 
model we have spherical symmetry and homogeneous and isotropic scattering conditions, therefore 
the same number of photons scatter out of as into the beam, and we can use the weight given in 
equation (4). In the case of infall we checked that equation (4) is still a good approximation in the 
regime in which we are interested. 

We evaluated the emergent intensity change at zero impact parameter (i.e. through the center 
of the gas sphere), where r = tq. In Figure 6 we show the intensity change A I for a static plasma 
for two optical depths and five temperatures. Non-zero optical depth causes only slight changes 
in the emerging radiation. Figure 7 shows the intensity change for a plasma with infall for two 
infall velocities and three plausible cluster temperatures. Only small changes in the spectrum are 
apparent, even with such large velocities. 
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One measure of the spectral deformation that has been used to quantify relative correction, 
and which is of use in determining the frequency at which to search for the KSZ effect, is the 
cross-over frequency. Challinor and Lasenby (1998a) and Birkinshaw (1998) suggested a linear 
expression for how the frequency changes with temperature, as 

X l 0 in = 3.830(1 + 1. 13 6 e ), (29) 

while Itoh et al. 1998 suggested a quadratic approximation 

Xg = 3.830(1 + 1.1674 0 e - 0.8533 0*) (30) 

for most cluster temperatures. In Figure 8 we compare our Monte Carlo results with these 
and other expressions that include relativistic corrections. Our Monte Carlo results for single 
scattering are close to those obtained by numerical integration of the collision integral (Nozawa et 
al. 1998). Relativistic corrections of third and fifth order (Challinor and Lasenby 1998a; Nozawa 
et al. 1998), or the linear approximation (Challinor and Lasenby 1998a; Birkinshaw 1998) are 
of varying accuracy in describing the curve: the linear and third order expressions give the best 
results, but extending the series to the fifth order is much poorer. This is a consequence of the 
asymptotic nature of the series, as emphasized by Challinor and Lasenby (1998a). 

Figures 9 and 10 show the cross-over frequency as a function of dimensionless temperature 
0 C for finite optical depth and infall velocity. Including a finite optical depth causes only a small 
change in the Xo(© e ) curve, and this change does not depend much on temperature. Based on 
our Monte Carlo simulations, we suggest the following approximation for the cross-over frequency 
for single scatterings in static spherical plasma for dimensionless temperature © e ^ 0.3: 


JX£(0 e ) = 3.827(1 + 1.2038 © e - 1.25670* +0.90980;!), (31) 

which fits better than 0.01 % in this range with a shift from the optical depth dependence 

A.Xq (to) = to (0.35 - 0.04416 ©e 0 ' 5 ), (32) 

which fits better than 0.005 % for 0.05 < 0 e < 0.3, and about 0.1 % for lower temperatures, and 
to < 0.5. 

Figure (10) shows results for our spherical models wiih gas infall. As for non- zero optical 
depth, additional energy transfers occur because of motion of the gas. As we would expect, at 
low temperatures the contribution to the electron velocity from bulk motion is comparable to that 
from thermal motion, and an enhanced frequency shift results, while at high temperatures this 
contribution becomes negligible. Based on our Monte Carlo simulations, we suggest the following 
approximation for the cross-over frequency shift in plasma with infall for 0.01 < 0 e < 0.3: 
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A^(0 e , (3 r ) = 0.224 ft 2 0-°- 7 . (33) 

This formula fits the cross over frequency v§ better than about half a percent. 

At small optical depth and bulk velocity we can assume that the shifts simply add, so that 
the final expression for the cross-over frequency becomes 


*o(0e, TO, ft) = X o *(0 e ) + AX 0 t (t 0 ) + A X$(Q e ,(3 r ). (34) 

In Figure 11 we show the cross-over frequency as a function of temperature and optical depth 
in the parameter range (0 e 0.04, tq £ 0.05) important for clusters of galaxies. From this 

figure we may come to the conclusion that, for clusters of galaxies, the optical depth effect on the 
cross-over frequency is more important than non-linear terms in the expansion in © e . We provide 
more accurate fitting formulae for this range. From our models we find 

Xg{& e ,T O ,0r) = 3.829(1 + 1.1624 0 e — 0.68948 0 2 ) +0.14 To + AXq (@ e , @ r ), (35) 

which fits better than 0.03 % for static models in this range of parameters. The optical depth 
independent first term is similar to the approximation provided by Itoh et al. (1998), which was 
obtained by numerically integrating the collision integral. 

The relativistic corrections to the kinematic SZ effect (cluster radial bulk velocity, v ra d) have 
also been found to be important (Nozawa et al. 1998, Sazonov and Sunyaev 1998b) The shift was 
found to be 


AX^ in (e e ,v rad ) = 300 


Vrad 


a l 

©„-©r 


+ 


02 


(©e ©mm)' 


+ 03 + 04 0 e + 05 © 2 


(36) 


where v ra d is the radial velocity, © min = 1.654 x 10 -3 , ai — 3.857 x 10 -3 , o 2 = -4.631 x 10 -6 , 
03 = 1.370 x 10 -2 , 04 = 1.014 x 10 -2 , and 05 = 0.01 (Nozawa et al. 1998), which agrees well with 
Sazonov and Sunyaev’s (1998b) result (expressed by a different fitting function). In Figure 11 , long 
dashed lines represent the results of this kinematic effect of cluster radial velocity of ± 100 kms -1 
(-Xo*(© e , 0 ) + AXo‘ n (0 e , Vrad)) Comparing our results for the shift in the cross-over frequency to 
results from relativistic kinematic effect, we conclude that a shift caused by an optical depth of 
0.01 is equivalent to a shift caused by a radial velocity of about v ra d = —10 kms -1 . 


We estimate the amplitudes of these effects on the hot cluster, Abell 2163, which was 
discussed by Holzapfel et al. (1997a). The maximum optical depth of the cluster is tq = 0.01, the 
temperature of the intracluster gas is close to 0 e = 0.03. Using our results for the static effect 
with the given temperature and maximum optical depth we get about 100 MHz shift to higher 
frequencies relative to the linear expression of Challinor and Lasenby (1998a). An infall velocity 
of /? r = 0.01 causes about an additional 15 MHz shift to higher frequencies relative to our result 
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for the static model. These shifts are small relative to the 20 GHz band width of the instrument 
of Holzapfel et al. (1997a) and the error in Hq and cluster radial peculiar velocity from ignoring 
their presence would be about 10 % (if the SZ effect is measured at x ^ 5) and 15 kms* 1 . By 
comparison, the component of primordial anisotropy in this scale corresponds to adding a velocity 
noise about ±200 kms -1 . 


4. Conclusions 

We investigated the effect of finite optical depth and bulk motion on inverse Compton 
scatterings in spherically symmetric uniform density mildly relativistic plasma. We assumed 
isotropic incoming radiation (CMBR), a relativistic Maxwellian distribution for the electron 
momenta, and scatterings in the Thomson limit. We demonstrated the usefulness of our Monte 
Carlo method for solving the radiative transfer problem, and calculated the static and kinematic 
SZ effects with different optical depth and gas infall velocities. 

The solution of the extended Kompaneets equation (with corrections up to the fifth order) is 
equivalent to a single-scattering approximation, and significant deviations from it occur for hot 
clusters and at high frequencies. These deviations may be as large as 5 % of the intensity change 
and neglecting them could cause about a 10 % error in the Hubble constant. A finite optical depth 
causes further small changes in the SZ effect: these changes may exceed the relativistic correction 
terms. For typical cluster temperatures, an accurate expression for the cross-over frequency as a 
function of temperature, optical depth, and bulk motion is (35). 

As it can seen from Figure 11, the cross-over frequency is sensitive to the cluster radial 
velocity, and less sensitive to the finite optical depth. Measurements of the cross-over frequency 
can, in principle, be used to determine the radial velocity of the cluster (e.g., as in Holzapfel 
et al. 1997b), with small extra corrections for optical depth and possible gas motion inside the 
cluster. However, the relatively strong variation of Xq with 0 e , compared to r or /3 r , suggests that 
the largest uncertainty will arise from the assumption of cluster isothermality, even if effects of 
confusion from primordial (and secondary) CMBR fluctuations can be excluded. 

Finally we note that our method can be extended to any geometry, density distribution and 
complicated bulk motion as desired, and may be used to study the SZ effect in high temperature 
plasmas with or without bulk motion. 
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Fig. 1. — The frequency redistribution function, P(s), for a static, single scattering case at 
dimensionless temperatures 0 e = 0.03 (a), or 0 e = 0.3 (b). The solid line shows the result 
from Rephaeli 1995b)’s semi-analytic method, the boxes with vertical error bars (too small to be 
visible) show results from our Monte Carlo simulations. 
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Fig- 2. — The intensity change AI/t (in units of i 0 — 2(k B T C B) Z /{he) 2 ) as a function of 
dimensionless frequency x = hi//(k B TcB) for dimensionless temperatures 0 e = 0.01, 0.02, and 0.03 
in static spherically symmetric models. The solid, dashed and dotted lines are results from Monte 
Carlo method (single scattering), from relativistic corrections to the Kompaneets approximation 
(Itoh et al. 1997) and the Kompaneets approximation (Kompaneets 1957) respectively. Note 
that at low temperatures, 0 e = 0.01 and 0.02, the single scattering Monte Carlo method and 
the Kompaneets approximation with relativistic corrections give very similar results (the solid and 
dashed lines overlap). 
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Fig. 3. — The frequency redistribution function, P(s), from Monte Carlo simulation of static plasma 
with uniform density at seven temperatures. The higher the temperature, the lower the peak and 
wider the function due to larger energy transfers from the hot electrons to the photons. 
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Fig. 4. — The inner part of the frequency redistribution function, P(s), for static clusters with 
dimensionless temperatures © e = 0.01, 0.02, and 0.03. The solid and dashed lines show the extremes 
of single scattering and scattering with maximum optical depth To = 0.1. More scatterings lead to 
more energetic photons, and hence more scattering from the line center to the high energy tail. 
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Fig. 5. — The single scattering frequency redistribution function, P(s), for infalling plasma with 
dimensionless temperatures 0 e = 0.01, 0.02 and 0.03. The solid, dashed and dashed dot lines show 
results for models with infall velocities /3 r = 0 (static), 0.03, and 0.05. 
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Fig. 6. — The intensity change AI/t (in units of io = 2(/cb7cb) 3 /(/ic) 2 )) as a function of 
dimensionless frequency x = hi// (hgTcB) for five dimensionless temperatures in static spherically 
symmetric models, for single scattering (solid lines) and tq =0.1 (dashed lines). 
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Fig. 7. The intensity change A I/r (in units of to = 2 (&eT(£#) 3 /(/ic) 2 )) as a function of 
dimensionless frequency x = hv/{kBTcB) for three dimensionless temperatures in spherically 
symmetric models with infall. Monte Carlo results are plotted using single scattering approximation 
for static gas (solid line) and gas with bulk motion (/3 r = 0.05, dashed line). 
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Fig. 8. — The cross-over frequency as a function of dimensionless temperature, 0 e . The solid and 
long dashed lines show our Monte Carlo results for single scattering, and the results of a numerical 
evaluation of the collision integral (Nozawa et al 1998). The dashed dot line represents a linear 
approximation (Challinor and Lasenby 1998a; Birkinshaw 1998). The short dashed and dashed dot 
dot dot lines use relativistic corrections of third and fifth order in © e . 




r 



Fig. 9.— The cross-over frequency as a function of dimensionless temperature, © e for static, 
spherically symmetric plasma. The solid and dashed lines are our results for single scattering and 
ro = 0.1, respectively. The dashed dot line shows a linear approximation (Challinor and Lasenby 
1998a; Birkinshaw 1998). Finite optical depths cause a shift due to multiple scattering which is 
almost independent of the temperature. 
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Fig. 10. — The cross-over frequency as a function of dimensionless temperature, 0 e , for our 
spherically symmetric model with infall. The solid and dashed lines are for infall velocities /3 r = 0 
(reference static model) and f3 r = 0.05. We used single scattering results to show only bulk 
motion effects. The dashed dot line shows a linear approximation (Challinor and Lasenby 1998a; 
Birkinshaw 1998). The effect of infall becomes inconsequential at high temperature since the infall 
speed becomes negligible relative to the thermal speed of the electron. 
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Fig. 11. — The cross-over frequency as a function of dimensionless temperature, © e , for several 
optical depths in a static spherical plasma. Our Monte Carlo results are from single scattering 
(solid line), To = 0.02 (short dashed line), 0.05 (dash, dot, dot, dot line), and 0.1 (dotted line). 
The dash-dot line, which is hardly distinguishable from the solid line, is the linear approximation 
of Challinor and Lasenby (1998a), and Birkinshaw (1998). As a comparison, we plot the (large) 
effect of a radial cluster velocity of ± 100 km s -1 (kinematic effect) with long dashed lines (Nozawa 
et al. 1998; Sazonov and Sunyaev 1998b). 





